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AUTOMATIC  ERR(Hi  CONTROL, 

THE  INITIAL  7ALUE  PROBLEM  IN 
ORDINARY  DIFFERENTIAL  EQUATIONS 


ABSTRACT 

This  report  describes  some  general  routines  which  will  automatically 
analyse  and  control  the  accuracy  of  numerical  procedures  on  high-speed 
digital  computers.  The  routines  described  are  applicable  to  any  system 
of  ordinary  differential  equations  for  which  a solution  is  being  sought 
to  satisiy  given  initial  conditions,  provided  that  the  numerical  method 
being  used  is  of  the  step-by-step  type.  Seven  such  numerical  methods 
are  listed,  along  with  the  resxilts  of  applying  the  automatic  control  to 
the  etud^  of  a spinning  top. 


IVeface 


mth  high  speed  digital  computers,  the  speed  of  running  of  a prob- 
lem is  naturally  much  less  critical  than  it  vss  with  slow  machines 
or  with  desk  computers.  In  fact  some  of  the  sinqpl<pst  computing  schemes, 
those  which  had  to  be  discarded  in  the  past  $s  beii^  unsuitable  for 
hand  computation  on  the  basis  of  time  consumed,  now  call  for  reconsid- 
eration. On  the  other  hand,  storage,  which  has  been  abundantly  avail- 
able for  band  con^uters  in  the  form  of  blank  paper,  published  tables, 
and  published  books  describing  computational  procedures,  becomes  a 
critical  and  controlling  consideration  for  high  speed  machines,  at  least 
at  present. 

However,  it  is  not  true  that  running  time  of  a problem  is  the  only 
use  of  time  that  needs  consideration.  One  cannot  ignore  the  time  needed 
in  coding,  code-checking,  error  analysis,  and  analysis  generally.  Now, 
since  the  high  speed  machines  can  make  multiple  choices  as  long  as  there 
is  a fixed  criterion  determining  such  a choice  at  each  encounter  of  an 
alternative,  any  sequence  of  arithmetic  operations  connected  in  such 
forms  as  "either  ....  or  ...  depending  on  >iiether  ...  or  not"  can  be 
coded.  The  machines  can  therefore  be  made  to  exercise  routine  judgements 
for  us.  A large  part  of  coding,  code-checking,  error  analysis,  and 
analysis  generally  can  be  so  formalized,  and  is  therefore  capable  of 
being  done  by  the  machines.  And,  indeed,  serious  attempts  are  currently 
being  made  to  get  machines  to  do  some  automatic  coding.  Fi|rthermore  a 
simultaneous  semi-automatic  coding  and  a saving  of  code-checking  time 
can  be  achieved  by  accumulating  a "library"  of  such  basic  subroutines  as 
polynomial  evaluation,  root  finding,  integrating  procedures,  etc.,  ary 
one  of  which  can  be  called  in  autanatically.  As  for  analysis,  one  could 
right  now  code  routines  performing  operations  on  power  series;  but  the 
bulk  of  "mechanical"  analysis  will  not  be  coded  until  machines  are  built 
with  a large  enough  memory  to  permit  us  freely  to  imitate  with  an  appro- 
priate use  of  machine  storage  that  extended  use  of  variables  in  ^ich 
they  are  called  "indeterminates"  (in  Qerman  - " Uribes tiimnte"). 

We  come  now  to  the  consideration  of  time  spent  in  error  analysis. 

No  extensive  computation  can  be  regarded  as  satisfactory  without  a sub- 
sidiary con5)utation  which  evaluates  its  accuracy.  In  planning  a coiq>U- 
tatlonal  schedule  an  important  use  of  such  error-analyses  is  to  deter- 
mine how  the  computation  should  proceed  (See,  for  example,  F.  J.  Ifcirray, 

CO)  if  the  accuracy  is  to  be  controlled.  We  are  Interested,  for  high 
speed  machines,  in  having  such  analyses  and  controls  coded  to  run  along 
with  the  basic  procedure  so  that  not  only  do  we  compute  our  error 
estimates  as  we  go  along,  but  we  also  save  the  time  that  would  be  lost 
in  stopping  the  basic  computation  to  perform  such  analyses  and  to  decide 
which  way  to  continue.  Ifcst  error  estimates  in  the  past  have  been 
deliberately  crude  in  order  to  permit  hand  computation  by  some  simple 
formula.  This  was  necessary  because  ,a  human  computer  was  justifiably 
reluctant  to  keep  track  of  the,  perhaps,  thousands  of  errors  accumulating 
due  to  approximating  formulae  and  methods  and  the  accumulation  of  round- 
offs. With  the  use  of  high  speed  machines  such  reluctance  should  disappear. 
It  is  conceivable  that  each  approximating  formula  can  be  coded  to  Include 


a parallel  computation  of  an  error  estimate  or  tiound  for  the  -erro^  due 
to  Initial  errors,  and'  to  the  so-called  "tnincation"  error  of  the'  formula. 
It  is  further  conceivable  that-^a  general  error  estimating  procedure  could 
be  coded  which  would  dq  three  things:  ^ 1 - it  wquld  assign  a separate 
mem0r7  position  to  the  error  of  each  confuted  qi^antity,  beginning  vdth 
those  designated  as  initial  errors,  2 - the  procedure  would,  for  each 
arithmetic  operation^  coiqjute  its  effect  on  the  errors  of  the  coiq}uted 
quantities,  including  an  pstimate  of  rounding-off  at  that  operation, 

3 - the  procedure  could  distinguish  among  certain  types  of  dlscrimiiw- 
ting  orders,  and  act  accordingly,  (Sach  discrimination  order,  according 
to  the  outcome  of  the  comparison  involved,  orders  either  a repetition 
of  a subroutine  or  the  beginning  of  another).  For  one  type  of  discrim- 
ination, error  estimates  are  to  be  accumulatpd;  for  other  types,  a prior 
analysis  must  be  made  available  for  proper  use  of  the  discrimination. 

This  analysis  provides  an  estimate,  with  accumulation  (see  S,  Qonn,  £23), 

Such  a procedure  could  be  coded  to  do  this  error  estimating,  in 
each  of  its  three  functions  mentioned  above,  either  by  finding  absolute 
bounds  or  by  finding  statistical  error  estiiiates  at  each  operation,. 

A general  error  estimating  procedure  of  such  a kind  pays  for  its 
generality  in  the  same  as  do  all  th»  so-called  "monitoring”  sub- 
routines, idiich  examine  orde^^s  before  carrying  t)iem  out,  namely,  by 
being  much  slower.  In  practice  it  would  seem  preferable  to  give  up 
this  complete  generality  in  favor  of  Separate  error  estimating  proce- 
dures for  large  classes  of  problems.  An  inverse  interpolation  method 
whereby  a machine  can  be  made  to  find  to  a high  accuracy  the  region  of 
indeterminacy  (due  to  roupd-off)  of  the  root  of  an  equation  ip  one 
numerical  variable  is  desoribed  in  ths  last  chapter  of  S.  Qorn,  Ci2» 

The  present  rfport  provides  such  an  eiroz*  estimating  procedure  for 
solving  the  ipltial  value  problem  with  systems  of  ordinary  differential 
equations. 


1.  Introduction 


We  will  now  discuss  a simple  scheme  for  estimating  and  controlling  ■ 
errors  in  the  solution  of  systems  of  ordinary  differential  equations* 

For  the  moment  let  us  restrict  ourselves  to  a single  equation! 

(l)  * f(t,  x). 

The  error  control  schemes  we  have  in  mind  assume  that  we  have  a 
single-step  numerical  integrating  procedure  so  modified  that  it  yields 
a reasonably  good  upper  and  lower  bound  for  x in  terms  of  the  value 
at  the  preceding  step.  We  will  give  several  examples  of  such  procedures 
below,  ^ter  we  will  see  the  effect  of  relaxing  the  bounding  conditions 
to  mere  error  estimating  conditions.  The  most  general  subroutine  dis- 
cussed below  applies  at  each  step  the  well-known  extrapolation  to  zero 
grid  method  (see  L.  F.  Richardson^  and  J.  A,  Qaunt,  CXD  to  obtain  such 
bounds..  It  in  turn  is  applicable  to  all  such  single-step  integrating 
procedures  as  those  of  Eulerand  Runge-Kutta. 


Figure  1 

Assuming,  then,  that  we  have  a number  of  such  procedures  for  single  step 
integration,  standard  modifying  procedures  can  be  made  to  apply  to  all 
of  them,  which  modifications  can  estimate  the  errors  sharply  at  each 
step  and  can  begin  the  operation  afresh  at  any  desired  point  with  modi- 
fied mesh  size. 

In  a region  containing  no  singular  points  (see  figure  l)  the  solu- 
tion curves  do  not  cross.  Hence  beginning  with  the  initial  condition 
x^  at  t^,  one  can  compute  the  bounds  x^^  and  3c^  at  t^j  then  beginning  with 

each  of  these  as  initial  conditions,  there  will  be  two  pair  of  bounds 
bracketing  the  corresponding  solutions,  so  that  the  maximum,  X2,  and 

minimum,  x^,  will  be  bounds  for  the  solution  sought.  This  process  can 

be  continued,  and  bounds  the  accumulated  error.  If  for  preassigned 
measure  of  accuracy,  C , the  z*egion  of  indeterminacy  of  x,  namely  x^^  - Xj^, 
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ever  becomes  bigger  than  the  procedure  can  be  set  to  divide  the 
mesh  slse  at  the  first  step  two  and  begin  the  Integration  over 
again.  This  will  be  continued  until  either  the  integration  is  conqpleted 
with  7 known  to  within  € or  else  the  mesh  site  has  gone  below  a given  6. 
(in  this  latter  case  the  integration  can  be  pernitted  to  proceed  with- 
out further  interruption).  The  total  tijne  taken,  with  all  these  nesh- 
halvings,  will  never  be  more  than  twice  the  time  required  by  the  last 
mesh  size  (namely,  that  mesh  deemed  necessary  accor^ng  to  the  error 
estimating  method  \ised)j  for  if  Tj  is  the  time  needed  for  the  last  run 

through,  the  preceding  run  through,  having  less  than  half  the  number  of 


mesh  points,  would  need  less  than 


1 


7 


4 

and  so  on  back. 


Consequently  the  total  time 

?•» 

as  stated. 


As  a result  of  this  time  consideration  all  that  is  needed  in 
choosing  an  initial  mesh  size  is  to  take  it  no  larger  than  that  re- 
quired to  be  printed  out,  and  preferably  no  larger  than  In  the  brac- 
keting conditions  discussed  below,  if  these  are  convenient. 


The  error  control' Just  described  can  be  coded  as  an  Independent 
procedure  applicable  to  any  Integrating  scheme  of  the  type  assumed. 
Its  logical  flow  chart  is  shown  in  Figiire  2,  beginning  at  input  i and 
ending  at  exit  e. 


If  an  initial  error  is  to  be  assumed,  we  can  even  choose 
2o  " *0  ^ at  t ■ tg,  where  ^ is  a bound  or  an  absolute  estimate  of 

the  initial  error » * * " 


Figaro  2 


Error  Control 
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At  the  place  marked  e in  figure  2 one  can  print  or  store,  as  one 
/wishes,  either  and  - 1/2  *vl 

^.i  ■ ‘^vjat  the  end  of  the  sequence,  or  the  complete  sequences, 

or  any  chosen  subsequences.  These  printing  options,  as  well  as  the  two 
boxes  directly  above  in  figure  2 are  already  included  in  the  general 
error  estimating  procedure  of  figure  11,  which  is  one  form  that  the 
section  from dc  to  p can  take.  Since  other  forms  are  possible,  the 
orders  from  Kto  S should  be  last  in  the  sequence  ^m  i to  e. 


2.  The  Initial  Value  Problem  on  High-Speed 
Digital  Computing  Machines 

2.1  - Formulation  of  the  Computational  Problem 

a numerical  solution  of  the  initial  value  problem  we  will  mean 
a set  of  numbers,  considered  to  be  approximations  to 

idiere  the  functions  x (t)  are  uniquely  defined  in  some  region  by  the 

XU  ^ 

order  system  of  ordinary  differential  equations  i 
dx  (t) 

(l)  •••»  n, 


and  the  .initial  conditions t 


(2) 


- X, 


1,  2^  n. 


V'  O'  V,0' 

k numerical  integrating  procedure  will  begin  at  t 


• X. 


V.o 


To  proceed  from  t 


t with 
o 


tj^  to  t ■ increments. 


A X^  will  be  computed  to  give 

j3)  X 


. ^ - X , -t-A  X ^ 

v,i+l  v,i  v,i 

as  approximations  to  x^(t)  at  t * 


Suppose  we  have  a numerical  Integrating  procedure  whose  "local 
truncation  error",  T*v  i»  ^ order  k ♦ 1,  In  other  words,  if  every 

numerical  operation  of  the  procedure  were  carried  out  with  infinite  pre- 
cision, it  would  yield  approximations,  X^  ^ such  that  if  j^(t)  is 

the  solution  of  equation  (l)  satisfying  the  initial  condition  x . (t. ) > 

XV^X  X 

^ then  at  t ■ t^^^^  we  should  have: 


(h) 

WJJ. 

/ 

idiere  h - t^^^  - Such  a numerical  integrating  procedure  is  said 
to  be  o4^  k*th  order. 


A ■ p (B)  will  have  the  standard  meaning,  nanely,  A/B  is  uniformly 
bounded  for  the  whole  range  of  variables  in  A and  B.  In  particular, 
for  our  application  A ■ 0 (h*)  means  kfV-  remains  boiuided  as  h-^0« 
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We  will 


show  in  section  2.3  that  equations 


inply 

(5)  (t^)  ♦ 0 (h^) 

in  a region  where  the  derivatives, 

are  uniformly  bounded. 


ST  ^ 


V ■ 


(U)  and  (2)  together 


1, 


2, 


• • e^  n 


When  a system  (l),  (2)  is  given,  one  is  faced  with  the  question 
of  how  a numerical  solution  in  the  above  sense  is  to  be  chosen  in  an 
efficient  way  to  attain  a desired  accuracy. 


We  begin,  in  the  next  section,  by  describing  a number  of  such 
numerical  procedures.  In  the  section  following  wo  discuss  the  growth 
of  errors  over  a number  of  steps.  In  chapter  3 we  will  describe  methods 
of  automatically  estimating  and  controlling  this  growth  of  errors* 


2.2  - Single-Step  Integrating  Procedures 

Particularly  convenient  for  machine  application  are  the  "single- 
step"'  methods.  In  these  the  computation  of  dX  . involves  no  values 

■#  ^ 

, of  the  X . prior  to  the  i 'th  step. 


The  convenience  of  such  methods  is  due  to  three  advantages  they 
possess* 


In  the  first  place  they  make  fewer  demands  on  the  memory  of  the 
machines  than  do  multiple-step  methods,  >rtiich  use  values  of  X . for 

j < i when  AX^  ^ is  being  computed.  Thus  larger  systems  of  equations 

can  be  accomodated. 


In  the  second  place,  and  even  more  important,  with  a single-step, 
method  the  Interval  of  integration,  h^  - t^^^  - t^,  may  bf  changed  from 

step  to  step  without  introducing  ary  change  in  the  formulae  for  computing 
A With  multiple-step  methods,  generally  complicated  interpolation 

procedures  are  required  for  a change  in  step-size.  An  appropriate  varia- 
tion in  the  step-size  aids  considerably  in  reducing  the  growth  of  error 
in  the  large,  as  will  be  remarked  in  Section  U below. 


This  definition  may  be  found  in  Bennett,  A.  A.,  Milne,  W.  E,,  and 
Bateman,  H.  - Numerical  Litegration  of  Differential  Equations  - 
Bull.  Nat.  Res.  Council,  0.  S,  92  (1933)  - pp.  51-8?. 


T 


n 


The  third  advantage  is  t,iat  single-step  methods  can  begin  directly 
from  the  initial  conditions,  whereas  multiple-step  procedures  require 
special  formulae  to  compute  the  first  few  values  of  X . , 

V/e  now  list,  both  for  reference  and  for  illustrative  purpoMS, 
some  of  the  best  known  single-step  methods.  They  are  all  applicable  to 
systems  of  first  order  ecjatit' jns.  For  brevity  we  will  use  the  vector 
notation: 

Ai  - V - 1,  2,  ...,  n|  - Agj^, 

for  quaiitities  with  the  subscript  v. 

Ml.  - "Euler's  method"  (see  levy  and  Baggott:  Q4.J  ) 


k » 1 - hj^fCtj^,  Ij^). 

M2.  - (See  Collatz  ) 

k - 2 h. 

^i  - ^ f(t^,  X^),  and 

Ali  - ^ 

M3«  - "Modified  method  of  Euler"  (Levy  and  Baggott,  loc.  cit.) 

k ■ 2 where  • • • ' 

- h.  . t^,  X^),  and 
‘ x '"h  * V ^i 

MU.  - "Kutta's  third  order  approximation"  (Levy  aiKi  Baggott,  loc.  cit.iy) 
k - 3 ^ (Aj^  ♦ 3 A^"),  where 

Ai  - h^f(t^,.X^), 

A^  - ^ i) > 

Ai  - h^f(t^  ^'i^‘ 

m 
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1 ■ 

/ 


M$.  - "Runge-Kutta"  (Collatz,  loc.  cit.) 


AXi 

■ S 

+ 2 A i 

t$  •)$ 

+ 2 A i + A i),  where 

- h^f(ti. 

x^). 

A'i 

- 

+ 7 V 

* 7 

- hjf(tj^ 

H 1 <). 

A*, 

* ^±i  Xj 

L )• 

m6.  - "Ruage-Kutta-Gill”  - This  is  a sli^t  modification  of  M5_which 
results  in  not  requiring  so  much  machine  storage,  (S.  Gill  - 

k - 4 4 i (^i  [2  - 

A . - h. f (t,,  I. ) 

• • ” h^f  (tj^  + h^j  ^ 

a"  . h^fct^  ‘ 5 V * 7 [-1  • /tl^i  * ' 

a"  . hjf  (tj^  {}  * y?|At)- 

M7.  - ''Kutta~N3rstr5m"  (B,  J.  h^strcm, 


k - 5 ^ (23A^  + 125  - 8l4i^^  * 125  A.p,  where 

■ h^f(t^,  X^), 

^'±  - h^f(ti  + 3 h^,  Xi  + 

2 6A'  +UA 

Ai  - h^f  (x-i  > I hi,  Xi  ♦ — 15 ) 
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where 


/2]a[) 


1 ■ 
/ 


i). 


Ai  - *■  h^i  Xi  ♦ 

*•  y 

^ iv  9 8 A'".  - 5b  A\  ♦ 9oA  . ♦ 6 A . 

Ai  - h,f(t,  . f h,,  X,  . — i i- i i ), 

V . 8A\  ♦ 10  A\  ♦ 36A  < ♦ 6 A4. 

A i - ^ 5 ‘'i*  * — h 


I 


i 

t 

I 


2,3  - The  Error  in  the  Large 

We  will  now  take  up  the  proof  of  2,1  (S)  in  the  case  where  all 
conpitations  are  assumed  to  be  carried  out  to  infinite  precision^  and 
will  then  see  idiat  modifications  are  necessary  due  to  the  fact  that 
they  are,  in  reality,  carried  out  to  a finite  number  of  places. 

As  in  section  2.1  we  let  ^ (t)  be  the  solution  of  the  system 
2,1  (1)  satisfying  the  initial  condition: 


(1) 


X .. 


Let  us  define  the  "error  functions"? 


(2)  (t)  - x^(t)  - x^^j,(t). 

Then  these  error  functions  satisfy  the  system  of  differential 
equations : 

%,i  ^ ^ 


(3) 


V' 

"ar 


dt 


■ * *n,i^* 

By  the  mean  value  theorem,  if  f are  of  class  C in  a region  R,  then 


(I:) 


- 5^1  Sr  I “n,i 


9f. 


where  the  partial  derivatives  (n  - l,2,...,n)  are 


evaluated 


^,i  “ ^*5  |i,v,i,-*> 


some 


f ■ 
/ 


(5) 


5 j.v.i  ■ ' 

1. 

Also  ^(t)  satisfy  the  initial  conditions t 

- \±- 


Nov  since  are  of  class  in  the  derlvatiTes  are  imifomly 

bounded  in  some  closed  subregion  S of  idilch  subregion  is  assumed  to 
contain  the  points  Supposoj  then> 

IsStI  ' ^ S. 

J 

It  follows  from  equation  (1*)  that 

(7)  |®v,i^*^|  - i /*n,i^*^j  * 

Let  us  define  E^(t)  bys 

(8)  E^(t)  - mM  * 

It  follows  from  equation  (?)  that 


(9) 


|bJ  ^ nLE^(t)  for  v - 1,2,,..,,  n. 


Now  clearly,  from  equation  (3),  since  f ciont  in' 

are  of  class  C for  t ^ t^t.^.,1  hence  E.  (t)  are  continuous  anSL^^e 
right-hand  derivatives*  D^S^(t)  are  piecewise  continuous  in  the  smoe 
intez^al.  Furthermore,  from  the  defining  equations  (8)  ve  have,  for 


each  t and  some  v varying  with  t In  t^i^  t»  t^^^^ 


that 


Hence: 
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»j^(t)  - B^(t) 
t - t 


D^Bi<t) 


lim 

^’-►t  ♦ o 


(10) 


I njs^Ml  $ -p  I I , ts 

It  therefore  follows  from, equation  (9)  that 

(U)  ' | D,Si(t)J  S'  nl*i(t),  ts 

It  is  now  possible  to  apply  the  following: 

Lenina:  (B*  Kajnke  -OO^  fflLlfssatz  2,  p.  93)  If  u(x)  is  continuous  on 

the  closed  interval  [^Sf  bj  and  has  a right  hand  derivative  for  a^*  b, 
and  if  there  exist  constants j M^O,  N5?0  such  that 

I K^u(x)  I ^ M I u(x)  j ♦ N,  then  for  any  two  nriaibers  x,  ^ in 
[a,  bJ  we  have 

|u(:c)|  ^ u (f)  -*l  * g (e«l^  -1). 

Taking  B^(t)  for  u(x),  nL  for  M,  and  0 for  N,  x « ^ “ "til 

and  h - t^^^^  - (in  the  remainder  of  this  section  h will  be  constant)^ 
we  have: 

(12)  ^ ®i<*i>®”“'* 

Now  from  equation  2.l(l+)  we  have  that  * 

^ *v,i^*i<M^  ° V (2)  and  (l), 

\i*i(  W ' ®»,i<  W " ° • 

Consequently,  by  (8),  (12),  and  (13): 

(lU)  *141^  Vl^^  S3^(ti)e*^  ♦ 0 (h*'*^). 


Applying  (Hi)  repeatedly  therefore  yields:  , 

(l5)  ♦ fl  ♦ e““*  ♦ e^*^n..oe^““i)o(h^^) 


But  since  I ^ • x^(t^)  Initially,  ^^(tQ)  ■ 0,  so  that  (l5)  becomes: 


(16) 


®i^l^^i+l^'^  ^nUi  T 
, ,e  -i  - 


n*om  the  fact  that  !♦!  » j ^^i+1  “ * 


hlh 


-1  ^ nUi  far  hUi  ^ o 


we  get: 
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'i 


\ 

I, 


(17) 


K. 


nDT 


-1 


0(h' 


k«l) 


whence 

(18)  E^(t^)  - m«  / x^(t^)  - I - 0(h^^ 

as  stated  in  equation  2.1  (5).  ^hls  conpletes  the  proof,  promised  in 
section  2.1,  that  if  a k^th  order  Integrating  procedure  Is  carried  oat 
with  infinite  precision  it  yields  an  approximation  in  the  large  of  order 
k in  a region  where  the  f have  unlfox^y  bounded  first  partial  deriv- 
atives . 


In  practice,  of  course,  conpitatlons  will  be  performed  with  some- 
thing less  than:  infinite  precision.  Real  numbers  will  be  replaced  b7 
rational  approximations  to  a fixed  nunber  of  places.  The  ftmctlons  f^ 

may  involve  infinite  series;  these  will  hitve  to  be  truncated  or  other- 
wise approximated  if  they  Involve  integrals,  they  may  be  approximated 
by  finite  sums;  if  th^  involve  roots  of  algebraic  equations,  these 
may  have  to  be  approximated.  In  general,  to  solve  digitally  requires 
reduction  of  all  computations  to  the  rational  operations,  and  the 
arithmetic  operations,  when  performed  on  two  s-place  numbers  do  not, 
in  general,  yield  exact  s-place  numbers;  the  result  must  be  truncated 
or  rounded  to  s-places.  Scaling  may  be  necessary  to  keep  all  nusber 
sizes  within  machine  range;  thls'may  result  in  the  loss  of  significant 
figures . * , 


To  see  how  this  finite  precision  affects  6\ir  accuracy,  lot  us 

represent  (for  the  remainder  of  this  section  only)  by  X . the  actual 

V j X 

approximation  correspohding  to  the  "infinite  precision"  approximation 

X . yielded  by  the  k»^th  order  method.  Similarly,  let  x . (t)  be 

y^x 

the  solution  of  equations  2.1  (l).  satisQrlng  the  initial  condition 
X .(t.)  ■ Then  in  place  of  2.1  (1*)  we  will  have 

*v,i^*i+l^  ” ■ T^v,i  * ^v,i. 


That  is,  the  local  error  will  consist  of  the  local  truncation  error 
TT  v,i  * 0(h  ) plus  an  error,  6 from  sources  mentioned  in  the 

preceding  paragraph.  For  example,  for  method  Ml  given  in  2.2,€^  ^ 
will  have  the  form  * 


(20)  ^ - h 

idiere  f^^^  ” ^2,±*****  ^n,i^  ■ X,  2,,..,  n, 

fy^j^  is  the  conQJuted  approximation  to  f ^ ^ h is  the  oomputed  approxima- 
tion to  h,  and . (^  is  the  rounded  moltiplication.  (Similar  expressions 


1 


IT 


i, 


f ■ 
! 


for^  . can  be  obtained  for  methods  M2  to  M7  in  2.2)*  In  this  case 

we  may  (in  the  manner  of  von  Neumann  and  Cteldstine,  decon^wse  ^ 
in  the  following  wayt  * 

* ^ Si  ■ “ ''.i  * ”*-i  * 

where  p . is  the  round-off  error  in  the  machine  product  of  S by  7 . . 

'^Vji  ! ^ ^ v,i 

For  machine  multiplication  with  radix  ^ and  rpundipg  to  s places, 

1 P 4 ^ !n  apy  case,  if  ^ .is  small  compared  to  the  local 

I ,v,i  > £ v,i 


trikication  error,  ^ .,  for  a given  h,  then  estimates  for  the  error 

X .(tj .,)  - I . -T  will  serve  as  well  for  estimates  of  x 4tt4j,)  - 
.,v,i  i+l  w,i+l  v,i  i'H'. 

However,  we  will  analyze  the  effect  of  t|ie  presence  of 

in  (19)  on  the  error  in  the  large. 

We  let  E .(t)  • X (t)  - x .(t).  E„  ^(t)  will  satisfy  the  initial 
condition 

(22)  - x,(V  - 

Corresponding  to  (3),  we  have 

(23)  *l»i#  ****  ^,i^* 

Wb  may  follow  through  (3)  to  (12)  exactly  as  before,  but  with  j^(t)  . 
and  x^,  ^(t)  Instead  of  E^  ^(t)  and  x^  ^^(t).  Wo  have  therefore, 

(2ll)  »i(W^*i^*i^ 

where  h - t^^^  - t^^  and  Bj^(t)  ■ muj  j^(t)  | . Nqw  from  (19)  wo  have 
\^*i^l^  ■ ^w,i*l  " %^^iel^  " S,i^*ii-1^  * ^v,i 

LotC  - max^  | € | ,T~  • «»en,  using  (2k), 

(26)  Vi(ti^lXii(ti)e'^*r*^‘. 


T 

I 


Now  initially  ^ where  is  the  error  li^  the 

initial  data  X^  T^eT^^'s  may  be  merely  the  roUnd-off  errors  com- 
mitted in  replacing  the  exact  initial  data  by  S-place  machine  nonberSf 
or  they  may  be  due  to  other  sources  of  error  such  as  measurement  in 
experiments  in  those  cases  where  the  initial  value  problem  directly 
represents  a physical  phenomenon.  In  any  event;  we  have 

(27)  (t^)  - - i,,0 

Letting  m^  l^yl  > repeated  application  of  (26)  yields 

(28)  (im”**..... , 


So  that,  for  fixed  t.^,  ■ t ♦ (i+L)h: 

' l-l  0 

nl(ti_^l-  t ) 

JEh — - ( t ♦ 6 ) ♦ ; 


We  recall  that is  the  local  truncation  error,  ■ 0(h^*^)  for  a k'th 
order  integrating  method,  and  where,  for  ML  for  example. 


(30) 

as  in  (21),  If  the  initial  data  X ^ are  taken  to  be  exact  (i.e  A.  • O), 

V^O 

and  if  S is  small  compared  to  T , then  for  a fixed  h a bound  such  as  (17) 
for  might  also  be  good  for  B^(t^).  However,  si?>pose  we  carry  a 

fixed  number,  s,  of  places  in  a sequence  of  computations  of  X^  de- 
creasing h in  the  sequence.  Then  lim  V * 0,  whence  by  (l7)  IS"E.(tJ  - 0 

h-»o  h-^o  ^ ^ 

On  the  other  hand,  (30)  shows  that  G cannot  be  expected  to  approach  0, 

whence  (29)  shows  that  (t^)  might  conceivably  grow  indefinitely  with 

decreasing  h.  In  fact,  then,  one  should  set  a lower  limit  on  h by  letting 
the  bound  in  (29)  be  minimum,  i.e.  take  no  h smaller  than  that  making 

V*€ 

- -r—  minimum.  Further  decrease  in  h would  not  only  increase  the 


labor  of  computation,  but  would  produce  worse  results  as  well.  Now  for 

k 


a k'th  order  integrating  procedure 


IT  ♦ £ 


for  sufficiently  small  h,  where  C and  g.  (■ 


behaves  like  Ch"’  ♦ ^o/b 

say,  11m  £ ) are  independent 
h-^o 


of  h.  This  ftinction  has  a minimum  when  h - (■  , idilch  should 


‘1 


be  the  smallest  value  \ised.  Note  that,  unlike  the  usual  criterion  (that 
It  Is  useless  to  take  h so  small  that  the  local  truncation  is  less  tbaii 
the  local  round-off)  this  permits  smaller  values  of  h,  namely  dgvn  to 

t.  ■ Ch  ■ “C  • 


2.i4  Variation  of  Step  Size  from  Step  to  Step* 

For  simplicity  the  analysis  of  the  preceding  section  was  carried 
through  for  step  size,  h,  fixed  throughout  the  range  of  integration* 
However,  as  was  mentioned  in  2*2,  it  may  be  more  efficient  canputattoa^Uy 
to  vary  the  step  size  within  the  interval  for  a given  problem*  We  havi 
just  discussed  a lower  limit  for  h*  Sax  upper  limit  is  required  to  kedj^ 
the  local  truncation  error  small  enou^.  In  a given  problem  we  ilnay  gain 
markedly  in  efficiency  by  varying  the  step  size  between  these. 

As  an  example,  one  procedure  for  varying  h which  is  sometimes  useful 
is  the  following:  (we  restrict  the  examnles  to  first  order  systems  - the 

extensions  to  hl^er  order  being  obvious)  for 


g - f(t,x)  , x(t^)  - Xq, 


choose 


"i  ■ *1*1  - ^1  ■ 


1 » (tj,  1^1 


vhere  is  the  numerical  result  for  x(tj^)  and  6 is  fixed  to  ma^fe  h 1^ 

between  the  limits  discussed  above.  This  method  spaces  (t^^,  X^)  at 

approximately  equal  intervals  of  arc-length.  Where  the  slope  o^  the 
integral  curVe  is  large,  l.e.  when  x is  changing  rapidly,  small’ steps 
are  taken,  whereas  for  slowly  varying  x,  where  the  curve  is  flat,  larger 
steps  are  taken.  When  f varies  widely  considerable  efficiency  can  be 
gained  by  not  keeping  the  step  size  always  small*  Essentially  the 
same  results  as  with  (2)  are  obtained  by  ■ ' 


n.  - , 

1 ♦ f(t^,  Xi)| 

which  is  con^nltationally  simpler* 

Another  procedure  one  mi^t  use  is  that  of  taking 


This  would  take  account  not  only  of  changes  of  slope  of  the  integral 
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curves  but  also  of  changes  in  the  rate  of  separation  of  neighboring 
integral  curves;.  This  method  can  be  modified  to  avoid  con5)uting 


where  X*  ^ say  X*  ■ X^  -ir  6, 


3*1  A Bounding  Procedure  for  Buler's  Method 

VTe  begin  hy  describing  hov  Boiler's  method  may  be  modified  to 
yield  a bounding  procedure.  Let  the  subscript  o applied  to  f refer 
to  evaluation  at  (x  , y ),  From 


y^  - f (x,y) 


vs  find 


y’.f' 

y*  - t'  - * 2 * Vy  * '4  • 


Figure  3 

Taylor  expansions  of  y and  f (x,y)  are,  therefore 


f (x,y)  - f ♦ hf  ♦ kt  „ ♦ 4 [h^  ♦ 2hkf  ♦ f . i 

o xo  yo  « 3GCO  xyo  779)  * 


uhere 


h - X - X 


’ k ■ y - y^« 

ti  o 


The  Buler  ipproximation  at  x^  ■ x^  ♦ h^^  is 

h • 2o  ♦ Vo- 
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A value  such  that  y^  is  bracketed  by  it  and  y^  is,  for  small 
enough  h^  and  an  unexceptional  Cx^,y^)  obtained  (see  figure  3)  by 
draving  through  (x^,y^)  a lino  parallel  to  the  direction  at  y^)i 
The  condition  on  (x  , y ) for  the  existence  of  such  an  hL,  and  the 

O 0 X 

limitation  on  h^  will  be  discussed  below.  Using  (3)  and  (5)*  (note 
that  k^  ■ ^1^0^ 


7l  - Vq  * h^f (x^,y^) 

- y j/f  +h,f  +h,f  f ♦ f «»>2hff  ^f  1 | 

^0  1 L o 1 xo  1 o yo  2 L 1 o xyo  1 o yyoj  J 

- y ^h^f  ♦hjf'#  i hhf"-  f^f') 

•^0  1 0 1 O e l'*o  yo  0 


From  (5)  and  (6)  we  find 

2 3 

(7)  I (?i  * y^)  - y„  * ^ *—  • 

Thus,  on  the  one  hand 

(8)  7 (y  - y ) - f{  + ^^^0^"^*  ' 


while, on  the  other 


»)  I ?)  - y,-  if 


Consequently,  if  there  is  no  point  of  inflection  at  (x  ,y  ),  l.e. 

o o 


<^0, 


j (^  + y^  ) is  a third  order  approxijnation  to  , while  ^ is 

only  oT  second  order.  This  is  what  we  meant  above  when  we  said  that 
the  modification  of  the  toiler  procedure  is  relatively  coarse.  Notice, 

however,  that  * %)  is  a third  order  approximation  ■nhether  ^ and 

bracket  3^  or  not.  This  well-known  approximation  Is  the  result  of  one 

iteration  in  the  so-called  modified  Baler  integrating  method.  It  Is 
also  known  as  the  Heun  method. 
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Now  let  us  find  the  condition  on  whereby  this  bracketing  is 

assured  at  one  step*  Suppose  that  in  a region  R containing  all  the 
points  (x,y)  in  which  we  are  interested. 


Vi 


where  1 and  J run  through  all  Integral  values  for  which  the  corres- 
ponding partial  derivatives  are  of  Interest.  (This  notation  is  simi- 
lar to  that  in  Bleberbach  Dlfferentlal-GQ.elchungen  p.  $$).  In  the 
present  case 

0 1^  1 ^2 

(12)  J O^j  ^2 

From  (2)  and  (11 ) we  therefore  have  in  R 

flf'l  ^ 2Ui 


f”  ^ 6L^M 


Then,  if  at  we  take 
' 0 


h ^ 


lOL^M 

we  will  have  bracketed  by  y^^  and  y^.  To  prove’  this^  observe  that 

by  using  the  inequality  for  f | and  f*‘  from  (13 ) in  (3)  and  (6),  ve  have, 
in  addition  to 

^1  - yo  * ^1  ^o» 

2 h^ 

that  y^  ♦ h^f„  * * h^f^  * T K* 

^o  ^ Vo  4^0  ■ 

p* 

if,  then  f^  >P,  it  follows  frCm 'several  ^plications  of  (lU)  that 

- y,  * Vo-^yo  * Vo  * r 
yi“yo  * Vo  * 

< y,  * Vo  * ‘i  ^o  - yi 


2lt 


On  the  other  hand^  If  f ^4.  0>  the  same  procedure  yields 

0 

2 

“i<  =^0  ‘ * Vo  * ^ ♦ ¥ ^o  - 

<y„  ♦ \f„  - lih%^  <3^  q.o.d. 

If  condition  (lli)  is  not  fulfilled  at  a given  point  for  a given 
mesh  size  h^  , we  have  remarked  as  a result  of  equation  (9)  that 

7 ^1+  is  still  a third  order  approximation  to  72*  Ejy  the  same 
token  equation  (8)  shows  that  ^ ^2  “ yi)  still  of  the  second  order; 
the  error  estimate  is  still  good  even  though  bracketing  may  be  lost* 

/ ?husj,  i£^the  initial  point  is  near  a point  of  inflection,  the  esti- 

mates are  of  the  same  order  as  above,  but  do  not  yield  bracketing* 

Suppose,  now,  that  (lit)  has  held  for  a number  of  steps  beyond  the 
initial  point,  so  that  at  3^,  y^  and  ^ oracket  y^^.  Then  even  though 

(ll;)  fail,  so  toat  bracketing  of  the  solution  curves  with  initial  con- 
ditions or  y^  or  f^ould  fail  at  ^ 7^41 

bracketed  by 


Figure  h 

In  other  words  (see  figure  U),  the  number  Of  steps  might  have  been 
sufficient  to  yield  a separation  from  the  actual  y^  bj'-  the  bounds,  which 

separ^ion  is  large  enough  to  counteract,  the  local  loss ' of  bracketing* 

In  still  other  words,  the  procedure  raif^t  have  had  a good  enough  running 
start  to  get  by  the  danger  point  without  loss . of  bracketing*  In  Uno 
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7 

/ 


i 


with  this  idea  we  may  proceed  farther.  We  can  give  a condition  on  the 
separation,  k,  between  and  y^  such  that  at  least  one  of  the  t|ouxds 

at  the  next  step,  beginning  at  y^,  still  ];vwnds  y^ 


Figure  $ 


X 

1 


1 


on  the  same  side.  !Die  condition  is 


(1^) 


k>- 


l-hL-h^L 


» 

(0^  |) 


(The  routine  mi^t  be  modified,  if  desired,  to  incorporate  (l5)  as  a test. 
This  would  ensure  bracketing  threu^out).  The  proof  is  similar  tb  that 
of  (lU): 

Since  y^  « y^  ♦ k, 

yj  - hf  (x^,  y^)  - y^  ♦ k + h |f^  ♦ kfy| 

^ • 

■ y.  * “o  * " * 

where  f is  evaluated  at  x and  some  y ■ 7q  ♦ k Sg,  0^92^1.  Thus, 
since  x^  ■ x^  ♦ h, 

" ^o  3^1  ' 

• y.  * >■  * *>  f'o-*  * [“o  * (1  * “y>]  V} 

■ 7o  » “o  * (f*  * f„f*)  « 1C  [l  * fl  ♦ hfy-] 
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* ■W’  •W'  / / 

where  f and  f are  evaluated  at  some  x + hfi, , y + k fio*  0^  1, 

xy 

2 

0^  ©2  — !•  But  y^  “ y^  + hf^  + ^ f ^ , where  is  evaluated  at 


k > max 


jh  f 

1 [*'  - 

2 (f*  ♦ f f*)]  ) 

X 0 y'J  ( 

» 

.« 

f T 1 

.1  + h f 

1 ♦ h f 

fl  + hf  1 

. y 

y 

L yJ  J 

(as  long  as  hL<-^“  ■-■) . By  using  (ll)  and  (l3)  on  this  we  get  (l5). 

q.e.d. 


Clearly  the  same  analysis  would  apply  to  y using  -k  instead  of  k. 

• • ^ 

• Note  that  if  we  begin  at  a point  of  inflection,  we  could,  to  pre- 
serve bracketing,  use  k in  (l5)  to  make  our  start  from  y - k and  y ♦ k* 

o * 0 

The  logical  flow  chart  for  the  above  modification  of  Euler's  method 
is  the  following; 


Figure  6 

Logical  Block  Flow  Chart  for  a Bounding  Procedure  on  Euler's  Method. 

27 


f ■ 


Hare>  because  the  F la  variable  in  application  the  coding  pf^ 
to08hould  be  last  in  sequence. 

3.2  A General  Error  Estimating  Procedure 

We  now  approach  a general  method  of  modifying  «iy  single-step 
integrating  procedure  to  produce  simultaneously  a hl^er  order  proce- 
dure and  an  error-estimating  procedure.  It  will  be  approached,  in  three 
stages . First  there  will  be  the  already  mentioned  method  of  extra- 
polation to  zero  grid.  Then  this  method  will  be  reorganized  to  be 
appll^  locally  at  each  step.  Finally  this  local  extrapolation  to 
zero  grid  will  be  modified  to  yield  bounds  at  each  step  if  the  grid 
size  is  small  enoughj  we  have  seen  in  the  introduction  how  this  would 
yield  error  bounds  at  each  step.  In  practice  it  will  not  matter 
-whether  the  grid  size  is  small  enou^  to  yield  strict  bounds,  since 
the  effect  will  be  to  produce  a higher  order  procedure  in  any  event, 
and  to  change  the  error  bounds  to  error  estimates. 

Let  uSj  then,  begin  by  recalling  the  method  of  extrapolation  to 
zero  grid,  with  a quick  proof  which  will  be  made  more  precise  in  our 
application. 


Figure  7 ‘ 

Extrapolation  to  Zero  Grid  in  the  Large 

Given  a single  step  integrating  procedure  idiich  ^)proximates  the 
exact  solution  (we  will  omit,  for  brevity,  the  subscrip>ts  for  the 
different  dependent  variables  until  we  set  up  the  finiG.  procedure) 
with  Initial  condition  z - Z_  at  t * t , x(t;  t , Z ) by  the  values 

0 O 0/0 


\ 


4 


(1) 


at  + h,  at  tg  * t^^  ♦ h,  etc.,  the  solution  process; 

was  run  through  again  with  grid  size,  say.  h/2  to  yield  the  valu,es  X^ 
at  t^  + h/2,  X^^^  at  t^,  X*  at  h/2,  Xy  at  tg,  etc.  (See  fig.  7). 

A better  approximation  to  x (t^j  t^,  X^)  is  then  given  by  the  formula 

2kj(2)  _ j(l) 

■ ..  — . We  have  seen  in  section  2.3  (equation  (l8))  that  a 

.2^-1 

”k’th  order  integrating  procedure",  irtiose  error  in  the  small  is  of 
order  k + 1,  actually  has  an  error  in  the  large  of  order  k.  Consequently, 
with  an  extra  differentiability  condition,  which  will  be  given  precisely 
below,  we  may  write : 


(1) 


x|^^-  X (t^j  t^,  X^)  ♦ 0 (h^) 


= X (t^;  t^,  X^)  ♦ C (t^,  X^)  h^  ♦ 0 (h^*^). 


Applying  this,  using  grid  size  h/2,  therefore  jdelds: 


(2) 


’'f  ^ ^ (‘i'-  ‘o'  V * “ 3 ° 


(2) 

It  follows  from  (l)  and  (2)  that  an  error  estimate  for  X^  'is  given  byt 


(3) 


k 

i . . C (t,,  .0 


A.  « —Tr 
1 2^-1 


while  a better  estimate  for  x (t.j  t , X^)  is  given  by 

gk  y(2)  yd) 

Xi  . if)  . Ai  - - X (t,,  t„,  I„)  . 0(h-*b. 


(u) 


“IT 

2*^-1 


i'  o'  o' 


This,  by  now,  classical  method  therefore  yields  either  an  error 

C2  ^ 

estimate  for  X,  ' or  a method  of  one  higher  order  for  x (t . ; t , X ), 
1 loo 


In  payment  for  this  one  has  two  separate  runs,  the  second  of  which  has 
twice  as  mariy  con5)utations;  thus  in  addition  to  the  con5>utations  (3) 
and  (I;),  the  running  time  is  multiplied  by  three.  Also  the  values  of 

the  first  run,  X^^^,  Xg^^,...  must  be  stored  until  such  time,  during  or 

after  the  second  run,  as  they  are  used  in  computations  (3)  and  (1|).  This 
last  storage  requirement  was  not  an  important  factor  in  hand  computations 
where  the  work  sheets  provided  the  storage.  In  high  speed  machine  com- 
putation it  is  a nuisance.  Let  us,  then,  as  a next  step  to  what  we  are 
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after  (an  error  bounding  or  estimating  method  to  <qppl7  to  the  resiiltlfcoiiig 
values)  see  vhat  the  effect  would  be  if  the  extrapolation  method 
reorganized  to  be  applied  at  each  step. 


Sxtrapolation  to  Zero  Grid  applied  Locally  at  each  Step 


In  this  procedure,  given  one  computes  and  froitHwha 

integrating  method,  then  one  computes  from  the  first  equation  in  0 (3)> 

after  which- one  computes  from  the  first  equation  in  (U),  and  onilss— 8 

ready  to  go  on  to  the  next  step  (see  Fig.  8).  Thus  the  whole  jmcocmi  *, 
from  “to  can  be  considered  a now  single  step  integrating  pTWh--** 

dure.  If, therefore,  we  show  that  locally  this  new  procedxire  is  bftisoae 
order  higher  than  the  one  it  is  applied  to,  the  general  result  in 
section  2.3  will  then  spply  to  yield  the  rest  of  equation  (U).  VeosB^t 
therefore  repeat  the  above  argument,  applying  it  locally,  and  we  tild  ~ 
this  opportunity  to  make  it  more  precise. 
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Figure  9 


We  assume  that  all  computed  values  lie  within  some  sin5)ly-connected 
region  containing  no  sin^lar  points  of  the  equations  2.1  (l),  and  further 

that  the  functions  f^,  in  these  equations  belong  to  class  C , i.e.  have 

bounded,  continuous  partial  derivatives  of  order  k+1  throughout  The 
first  k derivatives  are  needed  to  obtain  the  error  estimates  for  the 
k'th  order  integrating  procedure  so  that  we  may  write  equation  2,1  0*)|; 
the  extrd  derivative  permits  us  to  analyze  the  error  tarn  further  so  that 
we  have 


(5) 


k+1 


k+2. 


where  the  fundtions,  au?e  of  class  C^.  (It  is  now  essential  that  we 

exhibit  explicitly  the  subscript  v for  the  dependent  variables j see  Fig. 
9 for  notation,  otherwise) ; Applying  this  twice,  with  step  size  h/2, 
we  also  have: 


(6)  X 


v,i 


* 


h 


''r  ^li 


•'^ni^ 


hk*l 


o(h^^^). 


(7)  X 


(2) 

v,i+l 


h;  t. 


h Y# 


I*  ) 
•*'^ni^ 


Cv(ti 


h 


li 


, • « 


. jh 

ni-' 


k^l 

•3a 


0 (h^"^^)  . 
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Now,  definition. 


(8)  t^,  y^)  - y^, 

hence,  for  some  value  between  t^  and  ^ 7 

(9)  ♦ |)  V ^i*— ^ *ni^  • *vi  ^ I It  * ^i  * I ^ “ ^i  ^ 

Consequently,  applying  (?)  to  (6)* 

<,i  ■ V * oC**' 


whence 


(11)  ®v^^i  * 7*  ^j^i***>  ^n3[^  * ^i^**** 

(This  is  where  differentiability  assumptions  on  C^,  and  hence  on  the  f^ 
are  made). 

Applying  (ll)  to  (7)  now  yields  the  first  simplification  in 

(n\  # h^*^ 

(12)  - *v^*i  ♦ hi  *i  * 7*  *ni^  \i>**‘*  ^ni^p5T 

♦ 0(h*^*^). 

I 

In  order  to  be  able  to  use  this  with  (^)  in  the  same  wiqt  that  (2)  was  used 
with  (1),  we  must  analyze  the  first  term. 

Now,  in  general,  because  of  (8),  for  some  mean 't, 

(U)  x/tj  t„,  yj,...,  y„)  - y,  . 

Hence, 

• V * [ srI  t-ij  , 

Where  the  Krbnecker  6 means  1 if  v a |*  and  0 otherwise.  It  follows, 
using  (lit),  that 


x^(t 


GOC 

,(t;  t^,  y^,...,  y^)  - x^(t;  t^,  y^,..„  y^)  ♦ ^ ^ 

|S1  ^ 

Xy(ti  t^,  y^,..*,  y„)  ♦ (yy-yy)  ♦ (t-t^)  (y^^-y^^)  ^[sr]  tj 


H 

Now  let  t - t^  + h,  t “•  ■*'  7*  yj4  " 


and  7^  - V (tj  * §,-  V %>— > 

remembering  that  x^^t^+hj  t^  ♦ * 2* 

\(t^  ♦ ||  t^,  X^^)]  - x^(t^  •*-,h|  t^,  X^i>...,  by  definition, 


k*l 


and  that  ^ii*«**«  ■**  ® - oCh^"^^), 

and  t “ t^  ••  0 (h) . (l5 ) then  becomes : 

(16)  x^(t.  + hj  t.  + X^.,..„  X^^)  = x^(t.  + h;  t^,'x^.,..„  X^.) 


k+1 


^ni^  ^ * 0(h^*^). 


Substituting  this  into  (12)  now  yields 


*v^i+l  “ S^^i  V \i>‘**>  ^ni^  \i>***>  ^ni^"2k 


k+1 


♦ 0(h^*^). 


The  proof  now  proceeds  as  for  the  case  in  the  large:  from  (l7)  and  ($) 

(t) 

an  error  estimate  for  X^  is  given  by 

X^^^  - X^^^  k+1 


idiile  an  improved  approximation  for  + b|  t , X^^,  .o»,  is 


\,i+l  * ^v,i+l  ■^-^vji+l 


(19) 


„k^(2)  (1) 

v.1+1  ~ v.i+1 

2^-1 


V * 0(b^*h 


33> 


Equation  {h)  now  follows,.  In  the  large,  as  described  above..  In  this 
procedure  we  have  increased  the  order  of  the  given  integrating  procedure, 
but  are  now  without  any  sha^  error  estimates.  We  make  the  final  modifi- 
. cation  by  noting  that  locally  approximates  the  solution  to  one  order 

hi^er  than  ^PP^^o^cimates  the  last  correction;  thus,  for  a normal 

starting  point  we  should  have  ♦ | - X^^^  ♦ |AJ  and 

^ ^ ^ ♦ Aj_  *•  ( Axl  bounding  the  solution  above  and 

below  for  a sufficiently  small  h;  thenceforth  one  uses  X^  and  as 
initial  values  for  the  next  step,  computes  corresponding  and 

A i+i 

(see  figure  10).  In  other  words. 


(20) 

X * X -X 

- 

-v,o  v,0  V,0 

' - X^^' 

(21) 

A , 

2^-1 

__  X^2)  -"X^^;J 

(22) 

% 


(23)  X , - + A 4 - 1 A 4I  > 

— v,i  — v,i  I ' 

(2U)  T , ♦ A 4 ♦ 4I  r 

v,i  v,i  ^v,i  I t-i  v,i|  * 

are  the  formulae  for  the  new  procedure,  idiere  local  extrapolation  to 

zero  grid  on_X^  yield and  etc*  We  must  shoir 

that^  assuming  the  basic  integrating  procedure  and  the  aboYe  conqputatlons 
carried  out  with  infinite  precision,  we  will  have: 

(25)  o x^(t^i  t^,  x^)  ♦ o(h^*^), 

(26)  . 0(h‘'), 

* 

and  there  exists  an  h such  that 


(27)  x„(W  forh^h* 


Once  we  have  shown  these,  the  errors  6^  ^ described  in  paragrajii  2,3 

can  be  taken  into  account.  Equation  (27^  will  be  discussed  in  greater 
detail  below.  It  can  be  derived  immediately  from  (25)  and  (26)  as 
follows  J _ 

X . ♦ X . T - X 


Since  ■■-jA.,.g.7Zy.«.4,  . is  of  higher  order  than 


y,l  — v.i 


for  any  C > 0 there  exists  an  h such  that  for  h = t.  - - t^^h 

c j+1  j c 


(j  » 0,  1,...,  i-l)  we  have 


IX  . + 

1^ 


X . 


- *0’  ® 


Consequently  for  h ^ h^ : 


X 4 + X . 
v,i  — v,i 
1- 


T . - X , 

0 


+ c —**-5 — • 


Taking  c ■ 1 yields  (27). 


f ■' 

f 
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To  pravo  (2^)  and  (26)j  let  iis  apply  equation  (18)  to  the  defi- 
nitionSj  (2L)  and  (22): 

k^l 

By  applying  this  aid  il9)  to  the  definitions  (23)  and  (2U)  we  find: 

>^v,i^l  " *^*i+l*  V^Li****'-^i^  “ l * 0(h^) 

(29)  _ _ _ k*l 

^v,i+l  " ^^\*V  \*  h±***'*  *i^l  * 0(h*'^^), 


Novj  Just  from  the  cruder  information  in  (29)  that  the  errors  are  0(h  )t 
it  follows  by  use  of  the  result  in  section  2.3.  Just  as  equation  (5)  of 
2.1  followed  from  (U)>  and  from  definition  (2o)  above  that 

— v,i#l  ^o"  ^0'***>  *no^  * ■ ?^v,i^l* 

Equation  (26)  is  an  immediate  consequence  of  this.  Since  we  now  know 
that  and  differ  by  a quantity  of  order  k,  the  second  equa- 

tion of  (29),  using  the  fact  that  are  of  class  C,  becomes 

Iv,i+i  • V ^i>’*->  * K ^*i>  I ^ ^ 0(h^*^), 

Together  with  the  first  equation  of  (29)  this  yields: 

(31)  ^«i»l  ^ 0(h^^) 


The  proof  at  the  beginning  of  paragraph  2.3,  liien  applied  to 

Vt?ti.T^)«X',(t,  I,)  ■ 


®"vl 


”^8?  arlp«5“, 


instead  of  the  . (t)  in  (2)  of  2,3  and 
(U)  of  2.3  then  yields: 


I 


"v.i .\<vv v;^<vv-^? 


■ V *vo^  * 0(h**^), 

irtiich  is  (2^)  as  desired. 

The  logical  flow  chart  for  this  final  procedure  is  the  following! 

•» 

® Variables  t , H,  T,  t^, 


j 7li->-7V. 


X.v«>  Xv* 


V,  n, 

4 pf 


Xvl-^Xvi 
V = U."-,nrl, 


Sf—^S 

Kx~^A1: 


+ — ♦Qp 

Xvx“>iSyL 

^^|~^Xyx 

V‘l, 


' 

Figure  11 

Logical  Block  Flow  Chart  for  General  Error  Estimating  Procedure 


Mote  that  the  section  fl*omP  to  6 includes  the  integrating  pro- 
cedure* TUs  section  shotild  therefore  be  placed  last  in  a coded 
sequence,  for  «e  sav  in  section  2.2  how  much  it  can  vary  in  conplezity. 
By  the  sane  token,  in  the  coding  of  such  general  integrating  procedures, 
the  evaluation  of  the  f^  should  be  placed  last. 


If  the  procedure  is  to  be  adapted  tot  take  into  consideration  tbs 
errors,  £ „ described  in  section  2.3,  one  should  coiiq>lete  the  analysis 

of  the  ^ ^ in  the  particular  Integrating  procedire  and  system  of 

equations  by  finding  a boundt  end,  in  place  of  ^23)  and 

(2U)»  luw 


to)  ^ -JEif  ♦ ^ ^ -/4  J 


4lso,  errors  in  the  initial  conditions 
they  are  yj  ^ and  are  bounded  by  H^,  by 


(35) 


may  be  taken  into  account,  if 
using 


(36) 


♦ H 


Due  to  the  fact  that  these  last  local  errors  often  vary  in  sign,  in 
many  practical  problems  the  procedure  (20)  to  (2U)  more  than  suffice  to 
cover  them. 

At  this  point,  for  each  Integrating  procedure  an  analysis  like  that 
leading  to  (II4)  and  (l5)  in  section  3.1  is  possible  as  a supplement  to 
equation  (27).  »uch  analyses  increase  in  ccmplexlty  with  the  order  of 
the  single  step  integrating  procedure.  For  instance,  for  Euler's  method 
the  complexity  of  the  derivation  of  these  "sufficient  conditions  on  h for 
bracketing”  is  not  much  greater  than  the  derivation  in  section  3.1>  Ve 
will  not  give  it  here,  but  will  reiterate  that  most  of  the  time  sufficiently 
accurate  error  estimates  rather  than  strict  error  bounds  will  be  satisfactory. 


3.3  Example  and  Conclusion 


The  error-control  procedure  of  figure  2 and  the  error  estimating 
procedure  of  figure  11,  applied  to  the  Bunge-Kutta-Gill  procedure,  M6  . 
of  section  2.2,  using  constant  grid  size  were  applied  to  the  following 
problem,  a ginning  top  study  being  made  by  H.  need,  of  the  Computing 
Laboratory : 

u^  - 

"o  “ 

^0  " 

7o  * VUr 

z ■ 0, 

0 

t = 0. 

0 

As:  a check  on  the  results  the  following  identities,  corresponding  re- 
spectively  to  the  conservation  of  energy,  the  constancy  of  the  dyhi^cdl 
Arm,  -ashd  the ’ckj^hseryatia^  Spin^^n  be  used: 

2 2 

E = u + w ♦ y/U  = 1> 

2 2 2 

(38)  A ■ X + y + z S 1, 

S = ux  + ^ y + wz  3 1, 


• 1 „ 
u - ^ z, 

* > - I X, 

(37)  3c  - jjz  - wy, 

y » wx  - uz, 

. 1 

z - uy  - u^x, 


The  Problem  was  run  on  Edyac . The  maximum  error,  permitted  by  the 
control  was  set  at  5.  x 10"®,  and  the  minimum  step  size  (lixich  was 

never  needed)  was  set  at  6 » l/32j  the  initial  step  size  was  taken  at 
At  « 2,  since  the  first  set  of  print-outs  was  required  for  this  inter- 
val, The  control  reset ^t  to  1 after  one  step,  reset  At  to  1/2  arfter 
one  step,  and  computed  with  this  At  up  to  t « 2.$  when  the  error 
estimate  in  z went  out  of  bouhds;  it  finally  computed  with  At  ■ l/Ij  smd 
was  permitted  to  continue  to  beyond  t =»  3l|.  Finally,  the  control  was 
removed,  and  the  inte<?ration  plus  the  automatic  error  estimation  was 

run  with  step  size  At  » ^ up  to  T - 200  with  a running  time  of  about 

three  hours  and  print-outs  of  t,  u,  w,  x,  y,  z,  e(u),  e(w),  eOc),  e(y), 
e(z)  for  t = 0 (2)  10;  10  (lO)  200.  The  error  estimates  fluctuated, 
never  exceeding  (in  the  print-outs)  3$  x 10"7,  The  computing  time  was 
found  to  be  (not  counting  print-out  time)  9 mins.  30  secs,  per  interval 
of  t « 10,  an  average  of  7 secs,  per  step  of  Runge-Kutta-Qill  and 
automatic  increase  of  order  and  error  estimation  fOr  this  fifth  order 
system. 

•*I  \ 


?■ 

/ 


Saople  rfisults  are  tabulated  below 


t-30.  t - 1/8 

t-eoo.  ^t 

u -.81157642 

-.811576W; 

+.98238136 

e(u)  .00000222 

,00000008 

.OOOOOOilO 

w ♦*56^30lill4 

♦.56530liO7 

♦.13H7951 

e(w)  *00000262 

.00000009 

,00000099 

X -.7U349391 

-.7li3ii9l*31 

+.967li67li8 

e(x)  *00000076 

.00000002 

.000000U8 

y ♦.08709987 

♦.08710019 

+.07087767 

e(y)  .OOOOOlliO 

.00000005 

.00000257 

s ♦.6630U61i6 

♦*6630li58l 

♦.2li2861il2 

e(z)  .00000352 

.00000011 

.00000030 

E 1.00000002 

1.00000005 

1,00000610 

A 1.00000019 

.99999998 

.99999995 

S .99999995 

1.00000000 

1,00000023 

This  ezanple  tfould  indicate  that  vezy  little  practical  purpose  is 
served  in  deriving  formulae  corresponding  to  (Hi)  of  section  3*1  which 
would  provide  sufficient  coi^ditions  fnr  h to  be  small  enough  to  make 
the  error  estimates  actually  error  bounds.  However,  it  is  certainly 
conceivable  that  problems  will  arise  for  which  they  are  necessary. 

The  extrapolation  to  zero  grid  method  has  been  successfully  applied 
to  the  initial  value  problem  for  hyperbolic  partial  differential  equations. 
Se^  for  example,  F»  CUppinger  & N.  Qerber  W>.  C.  Carter  & G.  I. 

Spencer  seems  reasonable  to  expect,  tl^n,  that  the  automatic 

error  estimation  and  control  procedure  last  discussed  can  be  extended  to 
this  type  of  problem. 

4cknowledgements  are  due  to  Professor  A.  A.  Bennett  for  detailed 
editorial  mid  critical  advice,  and  to  H.  Bomanelli  for  the  coding  of 
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Bdvac  and  their  application  to  the  spinning  top  study. 
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